Skyrme Crystal in bilayer and multilayer graphene 
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The ground state of the two-dimensional electron systems in Bemal bilayer and ABC-stacked multilayer 
graphenes in the presence of a strong magnetic field is investigated with the Hartree-Fock approximation. Phase 
diagrams of the systems are obtained, focusing on charge density wave states including states with vortices of 
valley pseudospins (called a Skyrme crystal). The single-electron states in these stacked graphenes are given by 
two-component wave functions. That of the first excited Landau level has the same component as the lowest 
Landau level of the ordinary two-dimensional electrons. Because of this localized wave functions, the Skyrme 
crystal has low energy in this first excited level up to four layers of graphene, when the inter-layer distance is 
assumed to be infinitesimal. At the same time, bubble crystals are suppressed, so the phase diagram is diflrerent 
from that of a single-layer graphene. 

PACS numbers: 73.21. -b, 73.22.Gk, 73.22.Pr 
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I. INTRODUCTION 

Over the past decades, the conventional two-dimensional 
electron system (2DES) in semiconductor heterostructures in 
a strong magnetic field has been studied. It is found that 
electron-electron interaction brings various phases [e.g., frac- 
tional quantum Hall effect (FQHE) states''^ and charge den- 
sity wave (CDW) states-']. Similarly, graphene,"^ a flat sheet of 
carbons with a honeycomb lattice exfoliated from a graphite in 
2004,^ is under intensive investigation as a new 2DES whose 
electron has a valley degree of freedom {K,K') and a linear 
dispersion. Exhibition of FQHE^'^ and CDW-^ states in the 
new 2DES has been expected; FQHE was recently observed in 
a single-layer graphene (SLG).'"'" Furthermore, few-stacked 
graphenes have recently attracted attention due to their in- 
triguing properties: a band structure tunable by the number 
of layers and their stacking sequence, and a band gap control- 
lable by a perpendicular electrical field. '^ '-^ 

We focus on Bernal bilayer and ABC-stacked multilayer 
graphenes, which have chiral electrons concentrated on outer 
layers in low-energy states. In this paper, the CDW ground 
state of the 2DES in bilayer graphene (BLG) and multilayer 
graphene (M-LG, M > 3) in a strong magnetic field is stud- 
ied with the Hartree-Fock approximation. As candidates of 
the ground state, the electron (hole) Wigner crystal, electron 
(hole) bubble crystal, and Skyrme crystal are considered. A 
Skyrme crystal is the state with a topological texture of valley 
pseudospins. It has the lowest energy around a filling v = 1 in 
SLG.'' Our mean-field analysis predicts that the valley Skyrme 
crystal also has low energy in BLG and M-LG (M = 2, 3,4) 
at the first excited Landau level when the interlayer distance d 
is assumed to be infinitesimal. 

This paper is organized as follows. In the next section, we 
set up a low-energy effective Hamiltonian in the presence of a 
magnetic field for SLG, BLG, and M-LG. Then, the Hartree- 
Fock approximation is applied to the system of interacting 
electrons and self-consistent equations are derived. In Sec. |III1 
CDW states are introduced, including Skyrme crystals (more 
properly, the meron and meron pair crystal). As a preliminary 
calculation, the energy of a skyrmion (antiskyrmion) pair ex- 
citation at a V = 1 ferromagnetic state is evaluated and com- 



pared with that of a separated particle-hole excitation to find 
the condition where the Skyrme crystal has low energy. In 
Sec. |IV] numerical results are presented at the Landau level 
where the Skyrme crystal is expected. In Sec. [V] the validity 
of the results and a relation to experiments are discussed. 



II. MODEL 

A. Single-particle state in a magnetic field 

For a single-layer graphene (SLG), the low-energy effective 
Hamiltonian around the K valley is given by 



'Hf^ = Vf (p.,cr., + pyO-y:), 



(1) 



which has an electron with a linear dispersion.^-'"'"'^ Here, Vf 
is Fermi velocity, and cr,. and cr,, are Pauli matrices acting on 
sublattice (A, B) space. The Hamiltonian for the K' valley is 
'H^^ - -{'Hf'^y. In a perpendicular magnetic field B - 

Be,, using a magnetic ladder operator a = (tTx - i7Ty)lB/ 
{ti - p - eA), the single-particle Hamiltonian is written as 
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where Ib - ylWfeB is magnetic length. 
The eigenenergies are 



En = ±vf ^lefiBN, (A^ = 0, 1,...), 



(2) 
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where a positive (negative) sign is taken as the electron (hole) 
state. The coiTesponding eigenstates have the form 
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(A^ = 0), 





{<pN.x(r)) 

for the K valley. In the Landau gauge, (f>N,x(i') is given by 
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where L is the length of the system, Hn(x) is the Hermite 
polynomial, and X - kyfig is the guiding center coordinate. 
A macroscopic number of states with different X are degen- 
erated at one Landau level. A single-particle density on each 
sublattice Infig^'j^^^^x becomes broader as index increases. 
The states with broad density induce a bubble CDW and lose 
the profit to form skyrmions (see Sec. HiH i. 

Bernal-stacked bilayer graphene (BLG) and ABC-stacked 
multilayer graphene^" (M-LG) have low-energy quasiparticles 
with a dispersion E cc p'^, where M ( > 2) is the number of 
layers . '^I'^i^'i^^ In two adjacent layers, BLG and M-LG have 
interlayer hoppings between the sublattice /3 in the lower 
layer and the sublattice a in the upper layer In a magnetic 
field, the effective single-particle Hamiltonian for BLG and 
M-LG, which act on the outermost layer (top, bottom) space, 
have the same form 
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(M > 2), 



(6) 



where off-diagonal components represent M-times hopping 
via high-energy states in the inner layers. Its eigenenergy is 



En = +^wm y/N(N- 1) ■ ■ ■ (N - M + 1), 



(7) 



where Hcjm = f±( ^/invf/tJe)'^ ^ B'^'^. At zero energy, M 
Landau levels are degenerated. The corresponding eigenstate 
has the form 
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4'N-M,x(/) 

^ +4>N,xir) ) 

^ 
4>N.x(r), 



(N > M), 
(A^ < M), 



(8) 



Ib ~ 80 A, so we consider the vanishing Umit of interlayer 
distance as an approximate model. 

Apart from a constant kinetic term, the Hamiltonian for in- 
teracting electrons in M-LG (M > 1) is given by 



(T\ ,(T2,lT3,0"4 



(9) 



xy(|r-r'|)^^yr')^;yr), 
where the field operator i^JJ^^(r) is represented by 

X 

and Coulomb interaction is written as 
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er 



(10) 
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Here, e is the dielectric constant, cr = + represents valley K 
and K' , respectively, c'^^ is an annihilation operator of the 
electron at valley cr with guiding center X. A valley scattering 
term is relatively small when the cutoff q„,ax <k /ST is used j^'^-^^ 
since Ib » a { a - L42 A is the lattice constant), so we 
can ignore it as far as the CDW with magnetic length scale is 
concerned. Then, the Hamiltonian in q space is written as 

"^in' = ^2 Z Z Z Viq)p'^;i-q)p'^::iq), (12) 



N,X+igylfN,X-iq,Jl' 



(13) 



for the K valley. Notice that in the first excited state realized 
atN = M, the upper component of the wave function is given 
by <l>i),x{i'), which is the ground-state wave function of the or- 
dinary 2d electrons. 

At Landau level N < M, the valley degree of freedom co- 
incides with the layer degrees of freedom. If the layer degree 
of freedom in the M-LG model is regarded as the sublattice 
degree of freedom, Eq. (O is formally identical to that of the 
SLG model when M = L In the following, the case of M = 1 
(M = 2) represents the model of SLG (BLG). 



B. Hartree-Fock Hamiltonian 

Although the electron state is mainly considered in the fol- 
lowing, the hole state can be treated in the same way. In the 
present case of strong magnetic fields, electronic spins are 
completely polarized and the gap around the neutrality point is 
sufficiently large due to exchange enhancement; thus we take 
only one spin component and ignore the effect of Landau-level 
transitions. 

The interaction between the electrons in different layers is 
weakened by an interlayer distance d {~ 3.35 A for BLG). 
When B = 10 T, li is quite small compared to magnetic length 
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where Ln{x) is the Laguerre polynomial, and T^a-.e. is defined 
by A?+_| = A^_,j, = - M and A?+,x = A?_,| =" A^. When 
Landau-level index A^ is smaller than the number of layers M, 
the valley degree of freedom coincides with that of the layers 
(sublattices for SLG), and p|(r) = P\i^) - ™d P\^f) ^i^d 
p'^{r) are twice the amount of Eq. ( fT3T l. 

For the Hartree-Fock decoupling of "Kint, we assume the 
following order parameters: 

AT'iQ) - -jr n\ Ax ^ (15) 

X 

where Q is a reciprocal lattice vector of the CDW state. Then 
the Hartree-Fock Hamiltonian is given by^^-^Z 

= Z Z Z ^s:;^ (e)--'^'''A--(-e)c^Xc-;^_ 

Q X era-' 



Q X a-cr' 



N.X-' 



(16) 
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where Xj- - X + Qyf^/l. The Hartree-Fock potential consists 
of a direct term 



(17) 
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^"^.m-Y.^'MfiQ)^ (19) 
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Here, Jn{x) is a Bessel function of the first kind. 

The real-space density of electrons at valley cr and layer 
(sublattice) ^ is given by 

/'^W - ? Z (^TmLN^,iQ%me'Q'-'''''-l\ (21) 

The filling factor at Landau level is defined as 
y^ = ^e[0,2], 



(22) 



where - S /2nfg (S is a area of the system) is the degener- 
acy of Landau orbitals and A^^. is the number of electrons. The 
factor 2 comes from the valley degree of freedom. 



C. Green's function method 

To determine the order parameters self-consistently, the 
Green's function method is usedJi2i25iZ^ The single-particle 
Matsubara Green's function is defined by 



G^^^ {x,x',T) = -<rc-;,(T)4;,(0)). 



(23) 



The Fourier transformation 
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-(il2)QAX+r) 



(24) 



relates to the order parameters A?; °"(g) by 



(25) 



The equation of motion for G^'"^ (Q, r) is derived as 

Ji(icJn - H)G'^/' (Q, CJn) - Ti6Qfi6a-,cr' 

= ZZ^Mr;(e,e')Gf-'(e',a;„), (26) 

cr" Q' 

from the Heisenberg equation. Here, w„ is the Matsubara fre- 
quency for fermions, and the self-energy If^^j^iQ, Q') is given 
by 

cr" 

- x^ilQ - e'l)e-e^e'''/2A-' -(g - Q'). 

To solve this self-consistent equation, we diagonalize the self- 
energy matrix 



Q' 



_^.r/(e) 



(28) 



where (VJ, Vj) is the jth eigenvector with eigenvalue yj. 
The order parameters are obtained from the eigenvectors and 
eigenvalues 

Ar'(e) = Z ^(^^ - A')<(e)vr(o), (29) 

where f(x) is the Fermi-Dirac distribution function. The 
chemical potential fi is determined from 

^ A7(0) - ^ (0)VpO)f(yj -fi)^ VM. (30) 



Self-consistent equations are numerically calculated to yield 
the order parameters and the Hartree-Fock energy per particle 
for several CDW states introduced in Sec. |III1 

The order parameter sum rule at zero temperature^^ ex- 
tended to the case of valley degeneracy 



^^|A7'(e)|2 = A7(0) = y^ 



(31) 



is easily derived from Eq. ( |29] l. Here, is the contribution 
from valley <t electrons to the partial filling factor v^. This 
relation is used to check convergence of the results. 



III. CHARGE DENSITY WAVE STATES 

A. Valley skyrmion 

It is useful to map the valley degrees of freedom to a pseu- 
dospin.^ In this language, the components of the pseudospin 
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vector density P{Q) = P,{Q)x + Py(Q)y + PAQ)z are defined 
by 



PAQ) = 



PyiQ) = 



A+-(e) + A/(e) 



2i 



(32) 



(33) 



(34) 



In this paper, states with a topological pseudospin texture, 
which is called a skyrmion,^^ are considered. A skyrmion is 
a kind of spin texture usually used to describe a magnetic or- 
der and was first introduced in hadron physicsi^ For conven- 
tional and SLG 2DES, a state with aligned skyrmions, which 
is called a Skyrme crystal, has been shown theoretically to be 
the ground state around v = 1 



B. Preliminary considerations 

The excitation energy of a skyrmion in a ferromagnetic uni- 
form state at = 1 is evaluated to find the conditions on 
which the Skyrme crystal is preferred. When the interlayer 
distance c/ — » 0, the energy of a skyrmion (antiskyrmion) pair 
excitation As^ and the energy of a widely separated particle- 
hole pair excitation Ap^ are given by 

1 r°° 1 T 

^SK^-r q^V(q)[TM.N{q)f e-'^-'ll^dq, (35) 
4;r Jo 



1 r°° 

^PH^— qV{q)[TMAq)f e-^^'^l^dq, (36) 
27r Jo 



where the form factor 'FM,N(q) has the form 



(^) 



(A^ > M) 
(A^ < M) 



(37) 



The ratio of the two energies of the pair excitations is the same 
as that of single-particle excitations, from particle-hole sym- 
metry. For conventional, SLG, BLG, and M-LG (M = 3, 4, 5) 
2DESs, these energies at Landau level N < 5 are presented 
in Table H] It shows that skyrmion excitation is favored (1) at 
A^ = in all systems as single-particle wave functions being 
identical to conventional one, (2) at A^ = 1 , 2, 3 in SLG, and 
(3) at A^ = M in M-LG (M = 2,3,4). Thus we can expect 
that the Skyrme crystal becomes the ground state in these sit- 
uations around = 1. 

The reason why a skyrmion can be a low-energy excita- 
tion is the following.^"* The ground state at v^v = 1 is a pseu- 
dospin ferromagnetic liquid state. When a hole is introduced 
in this state without flipping pseudospin of other electrons, 
the charge density of the hole is given by an eigenstate of the 



TABLE I. Hartree-Fock quasiparticle and skyrmion (antiskyrmion) 
particle-hole excitation gaps (in unit of e- jelg ^njl) at Landau level 

and 



A' for 2DES in a conventional semiconductor structure (AJ,™'' 



A CO/IV 



and BLG and tri-LG (A; 



PH 



and hl-J^°), tetra-LG (A^„^° and l^f^), and penta-LG {h^p^ and 



-'SK 



^SK^^- ^^'^ conventional and SLG 2DES, the energies of the two 
excitations were compared by Yang et al^ The situation where the 
skyrmion is favored is emphasized by thick ">". 



N 



\conv. 
^PH 



Acanv. 



\SLa 
^PH 



\SLG 



vSLG 



A? 



BLG 



1 

0.75 
0.6406 
0.5742 
0.5279 
0.4927 



> 1/2 

< 0.875 

< 1.1328 

< 1.3418 

< 1.5522 

< 1.6834 



1 

0.6875 
0.5664 
0.5029 
0.4608 
0.4298 



> 1/2 

> 0.2188 

> 0.3301 

> 0.4097 

< 0.4754 

< 0.5328 



1 

0.75 
0.5977 
0.5029 
0.4528 
0.4187 



> 1/2 

< 0.875 

> 0.3770 

< 0.5151 

< 0.6181 

< 0.7048 



N M 



\3-LG 



^4-LG 



1 

0.75 
0.6406 
0.5498 
0.4660 
0.4223 



> 1/2 

< 0.875 

< 1.1328 

> 0.4448 

< 0.5808 

< 0.6829 



1 

0.75 
0.6406 
0.5742 
0.5285 
0.4406 



1/2 
0.875 
1.1328 
1.3418 
0.4870 
< 0.6284 



1 

0.75 
0.6406 
0.5742 
0.5279 
0.4963 



> 1/2 

< 0.875 

< 1.1328 

< 1.3418 

< 1.5522 

< 0.5390 



angular momentum, and is concentrated. On the other hand, 
when introduction of a hole is accompanied by pseudospin- 
flip of other electrons, many states with the same angular mo- 
mentum are connected by the Coulomb interaction, and the 
charge density of the hole has a wider distribution. This con- 
nected quantum state corresponds to a skyrmion. When the 
charge is confined like a wave function (/io,x given by Eq. (|5]), 
the formation of a skyrmion reduces the charge locality. At 
high Landau level A^, however, a wave function ^ai^ is intrin- 
sically broad, so the benefit to form skyrmions is lacking. In 
M-layered graphene, the spinor wave function has the local- 
ized component 0o,x at Landau level N - M, so it can drive 
the system to form skyrmions. 



C. Crystal structure 

When Zeeman energy for valley pseudospins does not ex- 
ist, a skyrmion splits into two merons (half-skyrmions). Four 
textures of a meron are possible from the two direction at cen- 
ter and the two vorticities. Charge density wave states where 
electrons, holes, or merons form a triangular or square lattice 
structure are considered. The order parameters A^'*^ (Q) are 
defined at points 



for a triangular lattice, (38) 



Q - (jGo, kQo) , for a square lattice, (39) 
where j and k are integers. The following states are assumed: 
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FIG. 1. (Color online) P,(r) for meron crystal at Vi = 0.86 in SLG. 



skyrmions per unit cell. Thus Qq is determined from 
the condition that the CDW has n (= 2) skyrmions in a 
unit cell: v/, = 2nf^n/s. The z component of pseudospin 
density in real space and the vorticity alternate from one 
site to the next (Figs. [T]and|2]l. The XY orientation of 
pseudospins has t/(l) symmetry (Fig.|2]i. The density 
distribution in real space is bipartite in layers. 

4. Meron pair crystal (MPC): a triangular lattice with four 
merons per unit cell. The merons are not equally spaced 
and bound into pairs. The go is determined from the 
same way as MC. The energy of MPC is similar to MC, 
and has slightly lower energy in the low quasiparticle 
density regime (close to = 1) in general. 



IV. RESULTS 
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FIG. 2. (Color online) The XY orientation of pseudospins for meron 
crystal at v\ - 0.86 in SLG. One meron crystal has two types of 
merons. 



\. Electron Wigner crystal (eWC) and «-electron bubble 
crystal (eBC«): a triangular or square lattice with one 
or n electrons per unit cell. The fundamental length in 
q space Qo is determined from the condition that the 
CDW has n electrons in a unit cell: v - InP^n/s. Here, v 
is the filling factor of electrons and s is the area of a unit 
cell. The low-energy state is pseudospin ferromagnetic 
because of the Pauli principle. 

2. Hole Wigner crystal (hWC) and «-hole bubble crystal 
(hBCn): a triangular or square lattice with one or n 
holes per unit cell. The fundamental length in q space 
Qo is determined from the condition that the CDW has 
n holes in a unit cell: v/, - InP^nj s. Here, v/, is the fill- 
ing factor of holes and s is the area of a unit cell. The 
low-energy state is pseudospin ferromagnetic because 
of the Pauli principle. 

3. Meron crystal (MC): a square lattice with four merons 
of charge -e/2 (v < 1) or e/2 (v > 1) per unit 
cell, equally spaced. A meron pair is equivalent to 
one skyrmion, so MC can be seen as a state with two 



We use the partial filling factor vm at Landau level A^; thus 
the total filling factor is given by v = AN - 2 + for any- 
layered graphene. Assume that Zeeman splitting is suffi- 
ciently large, so the phase diagram for e [0, 2] is identi- 
cal to that for vjv G [2,4]. Furthermore, the Hamiltonian has 
electron-hole symmetry around v^? = 1 , so the phase diagram 
for Vjv > 1 is caught by alternating particles for < 1 to an- 
tiparticles. Numerical calculation for MC and MPC are done 
at |1 - vai| > 0.06, since too many wave vectors are needed 
to get well-converged solutions ?A vn - 1 . In the following 
results, the energies with an accuracy of 10"*' are presented. 
Wigner and bubble crystals are calculated only in triangular 
symmetry; a square lattice generally has higher energy than a 
triangular one in a low quasiparticle density regime. 

It is difficult to get MC solutions close to - LO, 
so we extrapolate the order parameters of the MC solutions 
in vn < 0.94. To execute the extrapolation, the quantity 
'Fm,n{Q)^'^''^ (6) is fitted by a quadratic curve. The energies 
of the extrapolated MC states are represented in the following 
figures as a dotted line. It is noted that uncertainty remains in 
the extrapolation especially in the 4-LG case. It comes from 
the relatively low validity of fitting the order parameters which 
have an inflection point near - 0.92. 

The energies of the valley-concentrated hole Wigner crystal 
states which are quite accurately approximated by Gaussian 
form order parameters 



A-'-(O) = v^, ^'''''{Q ^ 0) = {vN - \)e 

A'^'*(e) = A^'"(e) = A*'^(e) = o. 



-e-/-/4 



(40) 



are represented as "GhWC" in the following figures. In the 
vanishing interlayer distance limit, the GhWC state with - 
(v^ - v^)/2 = +vai/2 is degenerated to the hWC solutions 
with < va,/2 (Figs, a Sill). 

In the following we show phase diagrams obtained by the 
present HF approximation in the whole range of v^. It should 
be remarked that the true phase diagram should contain re- 
gions of incompressible liquid states that cannot be obtained 
by the HF approximation. Thus the phase diagrams are partly 
incorrect. However, HF calculation gives qualitatively correct 
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results when the fractional quantum Hall states do not appear. 
This is established from comparisons with the results by the 
exact diagonalization method^^"^' or the density matrix renor- 
malization group method.^** In this paper we focus on the pos- 
sibility of meron crystals near = I, where the liquid states 
are not expected, but only the charge-ordered states compete. 
Therefore, the following discussion as to the realization of the 
meron crystal is reliable. 



A. Single-layer graphene 

The Hartree-Fock (HF) phase diagram of the 2DES in SLG 
has been obtained and compared with that of the conventional 
one.*^ It is shown that Skyrme crystals (MC and MPC) become 
the ground state around vai = 1 at Landau levels N = and 
1.'' Considering the excitation energies for SLG (TablelJl, it is 
also possible for the Skyrme crystal phase to occur at = 2 
and 3, but this has not been found yet in a mean-field calcu- 
lation. Although the same HF calculation had been done for 
SLG,^'^ we executed additional checks and investigated the 
higher filling regime. 

At Landau level N = 0, the phase diagram is the following: 
eWC for vo e [0.10, 0.50], hWC for vq e [0.50, 0.54], MC for 
vo e [0.54, 0.63], and MPC for vo e [0.63, 0.92].'^ 

At Landau level N = I, the phase diagram is the following: 
eWC for vi e [0.10, 0.50], hWC for v, e [0.50, 0.73], MC for 
vi e [0.73,0.84], and MPC for vi € [0.84, 0.92].'' The range 
of a skyrmionic (MC or MPC) phase is narrower than that of 
the N - case. 

At Landau level N - 2, the phase diagram is the follow- 
ing: eWCforv2 e [0.10,0.27], eBC2 for vz € [0.27,0.50], 
hBC2 for V2 € [0.50,0.73], and hWC for V2 e [0.73,0.94]. 
It is characteristic that 2-electron (hole) bubble crystals exist 
around V2 - 0.50. The skyrmionic ground state is not seen 
in the range V2 < 0.94. The extrapolating analysis, however, 
suggests that the hWC and MC state are almost degenerated 
in V2 e [0.94,1.0]. 

At Landau level N - 3, the phase diagram is the follow- 
ing: eWC for V3 € [0.10,0.20], eBC2 for V3 e [0.20,0.30], 
eBC3 for V3 e [0.30,0.50], hBC3 for V3 e [0.50,0.70], hBC2 
forv3 e [0.70,0.80], and hWC for V3 € [0.80,0.94]. The 
skyrmionic state does not appear in v < 0.94. 

Although the MC and MPC solutions are not found in < 
0.94 for SLG si N - 2 and 3, the Skyrme crystal is expected 
to have lower energy in the immediate vicinity of vai = 1 from 
the analysis in Sec. IIIIBI 
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for V2 e [0.50, 0.94]. The bubble state does not appear unlike 
the phase diagram at Landau level N = 2 for SLG. Figure |4] 
shows the energies in the area close to V2 = 1 . The extrapo- 
lated energy of the MC solutions have lower value than hWC 
forv2 € [0.94,1.0]. 



C. Landau level A' = 3 in trilayer graphene 



B. Landau level A' = 2 in bilayer graphene 

The HF calculation for degenerated zero-energy Landau 
levels = 0, 1 in BLG suggests that the Skyrme crystal states 
of real spin or orbital pseudospin occur ■'^ 

In what follows, the results for the first excited Landau level 
A^ = 2 in BLG are presented. Figure[3]shows the energies per 
electron for several crystal structures. It shows the following 
sequence of ground states: eWC for V2 e [0.10,0.50], hWC 



Figure |5] shows the energies per electron for several crystal 
structures at Landau level A^ = 3 in tri-LG. It shows the fol- 
lowing sequence of ground states: eWC for V3 e [0.10,0.46], 
eBC2 for V3 € [0.46,0.50], hBC2 for V3 e [0.50,0.54], hWC 
for V3 e [0.54,0.94]. Although the bubble states are found 
around V3 = 0.50, its range is narrower than that of the N - 3 
case in SLG. The skyrmionic ground state is not seen in the 
range v < 0.94. Figure|6]shows the energies near V3 - 1. The 
extrapolated states of the MC solutions have an energy close 
to that of hWCs in the vicinity of V3 = 1 . 
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FIG. 7. (Color online) Ground-state energy per particle (in units of 
e^/6/s) at Landau level yV = 4 in tetra-LG (cI/Ib = ). 
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FIG. 8. (Color online) Ground-state energy per particle (in units of 
e^/e/fl) around filling V4 = 1 at Landau level A' = 4 in tetra-LG 
{dllB = ). 



V. DISCUSSION 



D. Landau level A' = 4 for tetralayer graphene 



Figure I2] shows the energies per electron for various crystal 
structures at Landau level N - A for tetra-LG. It shows the fol- 
lowing sequence of ground states: eWC for V4 e [0.10,0.32], 
eBC2 for V4 e [0.32,0.50], hBC2 for V4 e [0.50,0.68], and 
hWC for V4 e [0.68, 0.94]. The bubble states are found in the 
broad range V4 e [0.32, 0.68]. The skyrmionic ground state is 
not seen in the range V4 < 0.94. Figure [8] shows the energies 
near V4 = 1 . The extrapolated states of the MC solutions have 
higher energy than that of hWCs in the vicinity of V4 = 1 . As 
previously mentioned, however, the extrapolation method is 
no longer valid in tetra-LG. According to the analysis in Sec. 
IIII Bl the Skyrme crystal is expected to have the lowest energy 
in the vicinity of V4 = 1 . 



The Hartree-Fock (HF) calculation suggests that the meron 
crystal (MC) phase appears at Landau level = 2 in BLG, 
and hole Wigner crystal (hWC) and MC states are degener- 
ated around v^? = 1 at Landau level A^ = 2 in SLG, A^ = 3 in 
tri-LG in the case of vanishing interlayer distance. Our calcu- 
lation strongly suggests that meron pair crystals (MFCs) will 
have lower energy in the immediate vicinity of vt^ - I for its 
triangular symmetry. 

Skyrme crystals (MC and MFC) have a charge distribu- 
tion and collective mode different from that of hWC, so these 
states can be distinguished by transport propertiesiSiii and a 
microwave absorption spectrum.^'*-' Furthermore, the CDWs 
exist in outer layers, so the local density of states (LDOS) 
can be measured in a spectroscopic manner.''^ In this paper, 
the unidirectional stripe phase is not considered. The stripe 
phase, however, also will appear around - 0.5 in BLG and 
M-LG, as SLG^ and conventional 2DES if fractional quan- 
tum Hall states are not realized. Such a state, if exists, will be 
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identified by anisotropic conduction4i 



Although we ignored the effects of disorder, a finite val- 
ley Zeeman energy, and Landau-level transitions, it is unclear 
how these affect degeneracy of hWC and MC around vn - I- 
In particular. Landau-level transitions in BLG and M-LG are 
larger than that of SLG under a magnetic field B ~ 10 T. 
The gap near a charge neutrality point is Ai - y/2hvf/lB ~ 
380 VWT K for SLG, A2 = y/lHcj,. ~ 45 x (B/T) K for 
BLG, A3 = y/enoji ~ 6.6 X (B/Tf'^ K for tri-LG, and 
A4 = yl7Ati(DA ~ LI X (BITf K for tetra-LG. The typical 
Coulomb energy is Ec = e^leh ~ 100 Vb/T K. For SLG, 
the ratio Ec /Ai - 0.39 is independent of field B. In this case 
it is shown that the Landau-level mixing does not change the 
CDW phase diagram (except the skyrmionic crystal)."*"* For 
BLG and M-LG, high magnetic fields are needed to achieve 
a comparable ratio: B ~ 70 T for BLG, 60 T for tri-LG, and 
50 T for tetra-LG. 



It has been pointed out that anisotropy in the pseudospin 
arises in the order of a/ls-^^'*^^^ In the multilayer graphene, 
finite layer separation also brings anisotropy. This anisotropy 
is quite small, since a < d <^ Ib, but may have some effect 
when the energies of two phases are quite close. We have 
done a calculation taking into account only the effect of d/ls 
as a preliminary investigation, and found that the finite d is 
slightly unfavorable for the Skyrme crystal. However, to ob- 
tain a definite conclusion for the effect of finite d, we need to 
take into account the effect of a also. Such calculation is left 
for future investigation. 
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